!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine model_variables (ro, mx, my, mz, bx, by, bz, en, gx, gy, po, x, y, nx, ny)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny

   real(8), dimension(nx), intent(in) :: x
   real(8), dimension(ny), intent(in) :: y
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, gx, gy, po

   real(8), dimension(nx, ny) :: pr
   real(8), dimension(0 : nx + 1) :: x_pre
   real(8), dimension(0 : ny + 1) :: y_pre
   real(8), dimension(0 : nx + 1, 0 : ny + 1) :: az_pre

   integer(4) :: i, j

   real(8) :: x0, y0, r2, r
   real(8) :: par_b, par_i, par_l1, par_l2, r0, h, dx, dy, b2, mag_e, beta, t_emerge_end, e_len, b_e

!----------------------------------------------------------------------------------------------------------------------------------|
!  magnetic flux function
!----------------------------------------------------------------------------------------------------------------------------------|
   dx = x(2) - x(1)
   dy = y(2) - y(1)
   x_pre(1 : nx) = x
   y_pre(1 : ny) = y
   x_pre(0) = x_pre(1) - dx
   x_pre(nx + 1) = x_pre(nx) + dx
   y_pre(0) = y_pre(1) - dy
   y_pre(ny + 1) = y_pre(ny) + dy

   t_emerge_end = 100.0d0
   e_len = 0.3d0
   b_e = 0.1d0

   do j = 0, ny + 1
   do i = 0, nx + 1
      az_pre(i, j) = 0.0d0
      if (y_pre(j) .le. 0.0d0 .and. abs (x_pre(i)) .le. e_len) then
         if (t .le. t_emerge_end) then
            az_pre(i, j) = b_e * cos (pi * x_pre(i) / e_len / 2.0d0) * t / t_emerge_end
         else
            az_pre(i, j) = b_e * cos (pi * x_pre(i) / e_len / 2.0d0)
         endif
      endif
   enddo
   enddo

   beta = 0.5d0
   h = 2.0d0
   r0 = 0.5d0
   par_b = 0.25628 / beta
   par_i = r0 / 2.0d0 / beta
   par_l1 = 1.0d0 / 2.0d0 / r0 / beta
   par_l2 = r0 / 2.0d0 / beta
   mag_e = 0.5d0

   do j = 1, ny
   do i = 1, nx

      bx(i, j) = 0.0d0
      by(i, j) = 0.0d0
      bz(i, j) = 0.0d0

!----------------------------------------------------------------------------------------------------------------------------------|
!     background field
!----------------------------------------------------------------------------------------------------------------------------------|
      x0 = -0.3d0
      y0 = -0.3d0
      r2 = (x(i) - x0) ** 2 + (y(j) - y0) ** 2
      bx(i, j) = bx(i, j) + par_b * (y(j) - y0) / r2
      by(i, j) = by(i, j) - par_b * (x(i) - x0) / r2

      x0 = 0.3d0
      y0 = -0.3d0
      r2 = (x(i) - x0) ** 2 + (y(j) - y0) ** 2
      bx(i, j) = bx(i, j) + par_b * (y(j) - y0) / r2
      by(i, j) = by(i, j) - par_b * (x(i) - x0) / r2

      x0 = -1.5d0
      y0 = -0.3d0
      r2 = (x(i) - x0) ** 2 + (y(j) - y0) ** 2
      bx(i, j) = bx(i, j) - par_b * (y(j) - y0) / r2
      by(i, j) = by(i, j) + par_b * (x(i) - x0) / r2

      x0 = 1.5d0
      y0 = -0.3d0
      r2 = (x(i) - x0) ** 2 + (y(j) - y0) ** 2
      bx(i, j) = bx(i, j) - par_b * (y(j) - y0) / r2
      by(i, j) = by(i, j) + par_b * (x(i) - x0) / r2

!----------------------------------------------------------------------------------------------------------------------------------|
!     image current
!----------------------------------------------------------------------------------------------------------------------------------|
      x0 = 0.0d0
      y0 = -h
      r2 = (x(i) - x0) ** 2 + (y(j) - y0) ** 2
      bx(i, j) = bx(i, j) - par_i * (y(j) - y0) / r2
      by(i, j) = by(i, j) + par_i * (x(i) - x0) / r2

!----------------------------------------------------------------------------------------------------------------------------------|
!     line current
!----------------------------------------------------------------------------------------------------------------------------------|
      x0 = 0.0d0
      y0 = h
      r2 = (x(i) - x0) ** 2 + (y(j) - y0) ** 2
      r = sqrt (r2)

!      if (r .le. r0) then
!         bx(i, j) = bx(i, j) + par_l1 * (y(j) - y0)
!         by(i, j) = by(i, j) - par_l1 * (x(i) - x0)
!         bz(i, j) = par_l1 * sqrt (2.0d0 * (r0 ** 2 - r2))
!      else
!         bx(i, j) = bx(i, j) + par_l2 * (y(j) - y0) / r2
!         by(i, j) = by(i, j) - par_l2 * (x(i) - x0) / r2
!      endif
!      ro(i, j) = 0.1d0

      if (r .le. r0) then
         bx(i, j) = bx(i, j) + par_l1 * (y(j) - y0)
         by(i, j) = by(i, j) - par_l1 * (x(i) - x0)
         ro(i, j) = 0.05d0 + par_l1 ** 2 * (r0 ** 2 - r2)
      else
         bx(i, j) = bx(i, j) + par_l2 * (y(j) - y0) / r2
         by(i, j) = by(i, j) - par_l2 * (x(i) - x0) / r2
         ro(i, j) = 0.05d0
      endif

      bx(i, j) = bx(i, j) + (az_pre(i, j + 1) - az_pre(i, j - 1)) / 2.0d0 / dy
      by(i, j) = by(i, j) + (az_pre(i - 1, j) - az_pre(i + 1, j)) / 2.0d0 / dx

!----------------------------------------------------------------------------------------------------------------------------------|
!     other variables
!----------------------------------------------------------------------------------------------------------------------------------|
      b2 = bx(i, j) ** 2 + by(i, j) ** 2 + bz(i, j) ** 2

      pr(i, j) = ro(i, j)

      mx(i, j) = 0.0d0
      my(i, j) = 0.0d0
      mz(i, j) = 0.0d0

      en(i, j) = b2 / 2.0d0 + pr(i, j) / gmm1

      gx(i, j) = 0.0d0
      gy(i, j) = 0.0d0

      po(i, j) = 0.0d0

   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
end subroutine model_variables
